arXiv: 1502.06742v 1 [math.OC] 24 Feb 2015 


VARIABLE DENSITY SAMPLING BASED ON PHYSICALLY PLAUSIBLE GRADIENT 
WAVEFORM. APPLICATION TO 3D MRI ANGIOGRAPHY. 


Nicolas Chaufferh 1,2 \ Pierre Weis s^\ Marianne Boucher^, Sebastien Meriaux^ and Philippe Ciuciu 11 


, 2 ) 


(I) CEA/DSV/I 2 BM NeuroSpin center, Bat. 145, F-91191 Gif-sur-Yvette, France 
(2) INRIA Saclay Ile-de-France, Parietal team, 91893 Orsay, France. 

< 3 > PRIMO Team, ITAV, USR 3505, Universite de Toulouse. 


ABSTRACT 

Performing k -space variable density sampling is a popular way of 
reducing scanning time in Magnetic Resonance Imaging (MRI). Un¬ 
fortunately, given a sampling trajectory, it is not clear how to traverse 
it using gradient waveforms. In this paper, we actually show that ex¬ 
isting methods min can yield large traversal time if the trajectory 
contains high curvature areas. Therefore, we consider here a new 
method for gradient waveform design which is based on the pro¬ 
jection of unrealistic initial trajectory onto the set of hardware con¬ 
straints. Next, we show on realistic simulations that this algorithm 
allows implementing variable density trajectories resulting from the 
piecewise linear solution of the Travelling Salesman Problem in a 
reasonable time. Finally, we demonstrate the application of this ap¬ 
proach to 2D MRI reconstruction and 3D angiography in the mouse 
brain. 

Index Terms — MRI, Compressive sensing, Variable density 
sampling, gradient waveform design, hardware constraints, angiog¬ 
raphy. 

1. INTRODUCTION 

Compressed Sensing (CS) provides a theoretical framework to jus¬ 
tify the downsampling of k -space (2D or 3D Fourier domain) in 
Magnetic Resonance Imaging (MRI). CS-MRI is usually based on 
independent random drawing of k -space locations according to a 
prescribed density. From recent theoretical works (3] 4|, one can 
derive an optimal sampling density i r that reduces at most the num¬ 
ber of samples collected in MRI without degrading the image quality 
at the reconstruction step BUD- In Q, simulations show that dis¬ 
tributions with radial decay (see Fig. [lja )) with full k -space center 
acquisition perform better in numerical experiments. 

However, such sampling schemes are not performed along con¬ 
tinuous lines and thus not physically plausible in MRI because of 
the constraints involved on the magnetic field gradient (magnitude 
and slew-rate). In (8), we have proposed a new approach to design 
continuous sampling trajectories based on the solution of Travelling 
Salesman Problem (TSP), as illustrated in Fig.fTJb). The specificity 
of this approach is that the empirical distribution of the trajectory can 
approximate any prescribed distribution 7r. Such a curve is called a 
7T-Variable Density Sampler (7T-VDS). Unfortunately, continuity of 
the sampling trajectory is not a sufficient condition in MRI and it is 
not clear how to design admissible gradient waveforms to traverse 
such a trajectory. 

To the best of our knowledge, the most efficient gradient de¬ 
sign strategies consist of finding an admissible parameterization of a 
given curve by using optimal control mm. However, the main draw- 
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Fig. 1. Example of 2D Variable Density Sampler, (a): 7r(k) oc 
1/1 /c | 2 as advocated in Q. (b): TSP-based sampling trajectory. 

back of this approach is that the traversal of high curvature parts of 
the curve is rather slow. In the extreme case of angular points such 
as in TSP-based trajectories, it leads to extremely large acquisition 
time. In (3, we have proposed a new method to design magnetic 
field gradients by projecting any parameterized curve onto the set of 
hardware constraints. This method allows one to change the curve 
support and thus to yield faster traversal time. 

The goal of this paper is then to prove that this projection algo¬ 
rithm enables to implement a TSP-based VDS on MRI scanners in 
reasonable scanning times, in contrast to optimal reparamerization. 
The first part is dedicated to summarizing the projection strategy in¬ 
troduced in (9). Next, we provide in Section [3] a proof-of-concept 
on retrospective CS simulations. In particular, we show that our al¬ 
gorithm yields faster sampling trajectories than the state-of-the-art 
for a given image reconstruction quality, and alternatively that if the 
scanning time is fixed, our method delivers improved reconstruc¬ 
tions. Finally, our strategy (TSP-based sampling + projection algo¬ 
rithm) is applied to 3D angiography of the mouse brain before and 
after contrast agent injection, to demonstrate its efficiency in terms 
of scanning time reduction while preserving the recovery of small 
structures such as blood vessels. 

2. GRADIENT WAVEFORM DESIGN USING A 
PROJECTION ALGORITHM 

In this section, we recall the hardware constraints as generally 
modelled in MRI and describe methods for designing gradient 
waveforms in order to traverse k -space sampling curves. 

The gradient waveform associated with a curve s is defined by 
g{t) — 7 _ 1 s(f), where 7 denotes the gyro-magnetic ratio |[K)j. 
The gradient waveform is obtained by energizing gradient coils (ar¬ 
rangements of wire) with electric currents. Owing to obvious phys¬ 
ical constraints, these electric currents have a bounded amplitude 



and cannot vary too rapidly (slew rate). Mathematically, these con¬ 
straints read: 

urn ^ G max and ||fl(t)|| ^ ^maxj Vt G [0, T\. 

A sampling trajectory s : [0, T] R d will be said admissible if it 
belongs to the set: 

S:={se (C * 2 ([0, T])) d , ||A(t)|| < a, p(t)|| < /?, Vi e [0 ,T]} 

In this paper, we limit ourselves to the so-called rotation-invariant 
constraints. Some hardware systems, where the coils are energized 
independently, enable considering rotation-variant constraints. In 
this case, the i 2 norm || • || is replaced by £00 norm. The differ¬ 
ences are discussed in m , but the two compared methods (optimal 
control and projection algorithm) are able to handle both kinds of 
constraints. 

2.1. State-of-the-art 

The question of finding jointly an accurate trajectory and the ad¬ 
missible gradient waveform to traverse it is a difficult issue that has 
received special attention in | Li 12]. The most classical approaches 
consist of fixing a curve support mm, or control points 1 101H31 . 
and finding an admissible parameterization afterwards. 

In particular, the most popular approach to design an admissible 
curve assumes the knowledge of a parameterized curve c : [0, T] — 
R d and consists of finding its optimal reparameterization by using 
optimal control mm. In other words, it amounts to finding a repa¬ 
rameterization p such that s = cop satisfies the physical constraints 
while minimizing the acquisition time. This problem can be cast as 
follows: 

T oc = minT' such that 3p : [0 ,T'] i-)> [0,T], cop g S. (1) 

The resulting solution s = c o p has the same support as c, which 
might be an important feature in some applications. The problems 
of this approach are: i) there is no control of the sampling density, 
especially in the high curvature parts of c where samples tend to 
agglutinate and ii) there is no control over the total sampling time 
T oc , which can be large if the trajectory contains singular points for 
instance (e.g., see Fig. |TJb)). These two drawbacks are illustrated 
in (9j. 

The next part is dedicated to introducing an alternative method 
relaxing the constraint of the curve support. 

2.2. Projection onto the set of constraints 

The idea introduced in (9) is to find the projection of the input curve 
c : [0, T] — R d onto the set of admissible curves S: 

s* := argminh|s - c||f. (2) 

where ||s — c||| := f^ Q ||s(£) — c(£)||J<it. For the sake of concise¬ 
ness, the theoretical grounds for the resolution and the key properties 
discussed below are given in (9). 

Resolution. Problem [2] consists of minimizing a convex smooth 
function over a convex set. In (9), we have proposed a fast itera¬ 
tive algorithm exploiting the structure of the dual formulation of 
which can be solved using proximal gradient methods m 
Key properties. The two main advantages of this method are that 
i) the sampling density of s* is close to the sampling density of c. 


In particular, if c is a 7T-VDS, the sampling density of s* is close to 
7i[j and ii) the acquisition time T is fixed and equal to that of the 
input curve c. In particular, the time to traverse a curve is in general 
shorter than with optimal reparameterization (T < T oc ). 

In the next part, we will emphasize that our algorithm enables 
to traverse VDS curves as depicted in Fig.[l|b) in a reasonable time, 
unlike optimal control-based reparemeterizations. 

3. CS-MRI SIMULATIONS 

In this part, we compare the time to traverse k -space along different 
trajectories using gradients computed either by the standard optimal 
control approach or by our proposed projection algorithm. For com¬ 
parison between sampling schemes, we work on retrospective CS, 
meaning that a full dataset has been acquired, and then a posteriori 
downsampling is performed. We compare the reconstruction results 
in terms of peak signal-to-noise ratio (PSNR) with respect to the ac¬ 
quisition time and to the “acceleration factor’^Jr. 

3.1. Experimental framework 

Data acquisition. The initial experimental setup aimed at observ¬ 
ing blood vessels of living mice using an intraveinous injection of an 
iron oxide-based contrast agent (Magnetovibrio Blakemorei MV1). 
Because of natural elimination, it is necessary to speed up acqui¬ 
sition to improve contrast and make easier post-processing such as 
angiography. The experiments have been performed on a 17.2T pre- 
clinical scanner which physical rotation-invariant constraints are, for 
all t G [0, T]\ 

\\g(t)\\ ^ 1 T.m -1 and ||g(t)|| ^ 5.3 T.m _1 .ms _1 . 

A FLASH sequence (Fast Low Angle SHot) has been used to reveal 
the T 2 contrast induced by the injection of the contrast agent (TE/TR 
= 8/680 ms). The sequence was repeated 12 times to improve the 
signal-to-noise ratio (SNR), leading to a total acquistion time of 30 
minutes to acquire the k -space slice by slice. The spatial resolution 
achieved is 90 x 90 x 180 pm 3 . 

Hypothesis. The aim of this paper is to prove that one can expect a 
large acquisition time reduction using partial k -space measurements. 
The time to traverse a sampling curve is computed satisfying the 
gradient constraints. To achieve a fair comparison, let us mention the 
additional hypothesis that our acquisitions are single-shot, meaning 
that the partial k -space is acquired after a single RF pulse. We did not 
take echo and repetition times into account to ensure the recovery of 
a TJLweigthed image. We only compare the time to traverse a curve 
using the gradients with their maximal intensity. We assume that 
there is no error on the /c-space sample locations. In practice we have 
to measure the three magnetic field gradients that are actually played 
out by the scanner to correct the trajectory and avoid distortions. We 
shall work on a discrete cartesian k- space, and consider that a sample 
is measured if the sampling trajectory crosses the corresponding cell 
of the k -space grid. Using this hypothesis, the estimated time to visit 
the 2D k -space is 110 ms. 

Strategy. We used the TPS-based sampling method (7) as input of 
our projection algorithm (see Fig.JTJb)), since it is a way of design¬ 
ing sampling trajectories that match any sampling density n. The 
latter is central in CS-MRI since it impacts the number of required 

Lhe closeness is quantified by Wasserstein transportation distance W 2 , 
see [9] for details. 

2 r quantifies the reduction of the number of measurements m. If the k- 

space is a grid of N pixels r := N/m is commonly used in CS-MRI. 



measurements muuia. To compare our projection method to ex¬ 
isting reparameterization, the proposed sampling strategy is: 

(i) Sample deterministically the k- space center as adviced in |[T4l l5l 
0, using an EPI sequence (see Fig.[2ja)). The scanning time can be 
estimated to 12 ms in 2D using optimal control. 

(ii) Select a density tt proportional to 1/| k\ 2 as mentioned in lH~5l 0. 

Draw independently points according to tt I and join them by the 
shortest path to form a 7T-VDS 0. 

(iii) Parameterize the TSP path at constant speed and project this 
parameterization onto the set of gradient constraints, or (iii bis) Pa¬ 
rameterize the TSP path using optimal control (the exact solution 
can be computed explicitely). 

(iv) Form the sampling curve, define a set Q of the selected sam¬ 
ples, mask the /c-space with D, and reconstruct an image using 
minimization of the constrained problem. Let F* denote the d- 
dimensional discrete Fourier transform and Fq the matrix composed 
of the lines corresponding to U. Denote also by an inverse d- 
dimensional wavelet transform (here a Symmlet transform). Then 
the reconstructed image is the solution of the problem: 

x* — Argmin (3) 

y= F h x 

An approximation of x* is computed using Douglas-Rachford algo¬ 
rithm 116|. Solving the penalized form associated with {3} might be 
addressed by competing algorithms (ADMM, 3MG); see H7) fora 
recent comparison. The reconstruction results could be improved by 
resorting to non-Cartesian reconstruction ED, which would avoid 
the approximation related to the projection onto the k- space grid. 


3.2.2. 3D angiography 

Using the same method as in 2D, namely TSP-sampling and pro¬ 
jection onto the set of constraints, we reconstructed volumes from 
3D k- space. In order to estimate the quality of the reconstructions, 
we compared the angiograms computed from the 3D images using 
Frangi filtering ED - The results are shown in Fig.|3|for acceleration 
factors r = 7.3 (Fig. [3jb,e)) and r = 17.4 (Fig. [3jc,f)) and com¬ 
pared to the angiogram computed from the whole data. 

Using the strategy described in Part |3.1| the time to traverse k- 
space would be 3.53 s (full acquisition), 3.15 s (r = 7.2) and 0.88 s 
(r — 17). The main drawback of TSP-based sampling schemes is 
that the time reduction is not directly proportional to r, in contrast 
to classical 2D downsampling and reading out along the third di¬ 
mension. Nevertheless, if the number of measurements is fixed, the 
TSP-based approach leads to more accurate reconstruction results 
since the sampling scheme may fit any density m. 

Angiograms shown in Fig. [5] illustrate that one can reduce the 
travel time in the k- space and still observe accurate micro vascular 
structure. If r — 7.3, time reduction is minor (about 10% less), 
but the computed angiogram is almost the same as the one obtained 
with a complete k- space. It is interesting to notice that with a higher 
acceleration factor (r = 17.4), the acquisition time is reduced by 
75%, but the computed angiogram remains of good quality. The 
angiogram appears a bit noisier, especially in the pre-injection set¬ 
ting (Fig.[3jc)), but the post-injection image allows recovering Willis 
polygon and most of the major vessels of the mouse brain (Fig.[3jf)). 


3.2. Results 

3.2.1. 2D reconstructions 

In this experiment, we considered a 2D k- space (d = 2) correspond¬ 
ing to an axial slice. We considered five sampling strategies, de¬ 
picted in Fig. [2|first row): a classical EPI coverage used as refer¬ 
ence (a); a TSP-based sampling trajectory parameterized using opti¬ 
mal control (b); two projected TSP-based trajectories, one with the 
same number of samples collected as in (b) (r = 11.2) (c) and the 
other with the same scanning time as in (b) (62 ms) (d); a variable 
density spiral trajectory for comparison purpose in terms of time and 
sampling ratio (e). 

As expected, the reconstruction results shown in Fig. |2[g,h) are 
really close, since the number of collected samples is the same, and 
the sampling densities are similar. However, in this comparison the 
gain in traversal time is significant (one half). In contrast, the longer 
and smoothed TSP depicted in Fig. [2|d) allows us to improve im¬ 
age reconstruction (1 dB gain) as illustrated by Fig. @i) while keep¬ 
ing the same acquisition time as in Fig.[2|b). For comparison pur¬ 
poses, we implemented spiral acquisition which consists of replacing 
steps (ii)-(iii) in the above mentioned sampling strategy by a spiral 
with density proportional to l/|/c| 2 , projected onto the set of con¬ 
straints. This strategy doubles the acquisition time (118 ms com¬ 
pared to 62 ms) whereas the acceleration factor was larger (r = 7.5 
vs. r = 6.6). In this experimental context (regridding and variable 
density spiral), the spiral is not appealing compared to EPI acquisi¬ 
tion, since it is time consuming and degrades the image quality. 

In each of these reconstructions, the major vessels can be recov¬ 
ered, although the smallest ones can only be seen for r < 8. Fi¬ 
nally, the best compromise between acquisition time and reconstruc¬ 
tion quality is achieved using the specific combination of TSP-based 
sampling and our projection algorithm onto the set of constraints 
shown in Fig.[2jd). 



Fig. 3. Angiograms computed from full /c-space pre-(a) and post- 
(d) injection data. Angiograms computed from pre-(resp., post-) 
injection data for decimated k- space with r = 7.3 (b) and r = 
17.4 (c) (resp., (e) and (f)). 










Fig. 2. Full /c-space acquisition with an EPI sequence (a) and corresponding reference image (f). Comparison between an exact parameter¬ 
ization of the TSP trajectory (b) and projection from TSP trajectory onto the set of constraints (c),(d). In experiments (b,c), the number of 
measured locations is fixed to 9% (r = 11.2), whereas in (b,d), the time to traverse the curve is fixed to 62 ms. (e): Spiral trajectory with full 
acquisition of the k- space center, (g-j): Reconstructed images corresponding to sampling strategies (b-e) by solving Eq. {3}. 


4. CONCLUSION AND FUTURE WORK 

In this paper, we have shown that the projection algorithm intro¬ 
duced in 0 is a promising technique to design gradient waveforms. 
In particular, it allows us to design gradient waveforms in order to 
implement TSP-based VDS, while existing gradient design methods 
lead to extremely large acquisition time. We are currently imple¬ 
menting these waveforms on actual scanners to validate our method 
on a real CS-MRI framework. 
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